function ptr (L)

   ## usage: ptr (L)
   ## 
   ## process predictors (analyses and GCM simulations)

   if nargin == 0, L = true ; endif

   global PAR PRJ GVAR ADIR SDIR CDIR isANA isCUR keof

   dclim = fullfile(CDIR, "clim") ;

   cfile = fullfile(CDIR, [strrep(PAR.ana, "/", "_") ".pc.ob"]) ;

   if !isCUR && !isnewer(cfile, PAR.mfile)
      wsim = PAR.sim ;
      init(["SIM = \"" PAR.sim{2} "\";"]) ;
      ptr(L) ;
      PAR.sim = wsim ; clear wsim ;
      init(["SIM = \"" PAR.sim{1} " " PAR.sim{2} "\";"]) ;
   endif

   ofile = fullfile(SDIR, [strrep(PAR.ana, "/", "_") ".pc.ob"]) ;
   if (exist(ofile, "file") && !L), return ; endif

   f = dbstack ; printf("\n\t###   %s   ###\n\n", f(1).name) ;

   PAR.prbwgt = [1 0] ; PAR.Pzero = 0.1 ;

   if exist(file=fullfile("data", PRJ, PAR.ptr, PAR.sim{1}, "ptr.ob"))

      printf("ptr: <-- %s\n", file) ; 
      load(file) ;

      if isempty(GVAR)
	 GVAR = [fieldnames(ptr) fieldnames(ptr)]' ;
      endif

   else

      ptr = read_ptr(fullfile("data", PAR.sim{1})) ;

      if isempty(GVAR)
	 GVAR = [fieldnames(ptr) fieldnames(ptr)]' ;
      else
	 printf("ptr: --> %s\n", file) ; 
	 save(file, "ptr") ;
	 ptr = SelFld(ptr, GVAR(2,:)) ;
      endif
      printf("ptr: --> %s\n", file) ; 
      save(file, "ptr") ;

   endif

   ## regrid and normalize
   if !isnewer(zfile=fullfile(SDIR, "ptrz.ob"), file, dclim) || length(dir(dclim)) == 2

      pkg load general nan optim statistics

      for iptr = 1:columns(GVAR)

	 k = GVAR{2,iptr} ; v = ptr.(k) ;

	 v = to_grid(v, k, iptr) ;

	 if PAR.ptrout
	    lon = v.lon ; lat = v.lat ;
	    save(pfile = fullfile(SDIR, [k ".geo"]), "-text", "lat", "lon") ;
	    csvwrite(pfile = fullfile(SDIR, [k ".csv"]), [v.id(:,1:3) v.x]) ;
	    printf("ptr: --> %s\n", pfile) ; 
	 endif

	 v.isP = !isempty(strfind(tolower(v.name), "precip")) ;

	 ## transform to N(0,1)
	 cfile = fullfile(dclim, [k ".ob"]) ;
	 v = do_norm(v, cfile) ;

	 ptr.(k) = v ;

      endfor
      save(zfile, "ptr") ;
      printf("ptr: --> %s\n", zfile) ; 

   else

      load(zfile)
      printf("ptr: <-- %s\n", zfile) ; 

   endif

   ## EOF reduction \\
   ptr = common_def(ptr) ;

   pkg load nan ;
   PC.id = PC.x = PC.J = [] ; wptr = {} ; keof = 0 ;
   for iptr = 1:columns(GVAR)

      if !isANA
	 printf("ptr: pairing %s with %s\n", GVAR{[2 1],iptr}) ;
      endif
      v = ptr.(GVAR{2,iptr}) ; wptr = {wptr{:}, GVAR{2,iptr}} ;

      [pc v] = do_eof(v, GVAR{:,iptr}) ;

      PC = upd_PC(PC, v, pc, iptr) ;

   endfor
   PC.ptr = wptr ; [PC.nr PC.nc] = size(PC.x) ;

   if isANA
      PCw = selper(PC, sdate(PC.id, PAR.cper(1,:), PAR.cper(2,:))) ;
      G = covadj(PCw.x) ;
      printf("ptr: --> %s\n", Gfile = fullfile(ADIR, "G.ob")) ; 
      save(Gfile, "G") ;
   elseif isCUR
      PCw = selper(PC, sdate(PC.id, PAR.cper(1,:), PAR.cper(2,:))) ;
      G = covadj(PCw.x) ;
      printf("ptr: --> %s\n", Gfile = fullfile(CDIR, "G.ob")) ; 
      save(Gfile, "G") ;
   endif

   if strfind(PAR.ptradj, "mlt")
      PC = mltadj(PC) ;
   endif 

   save(ofile, "PC") ;
   printf("ptr: --> %s\n", ofile) ; 
## EOF reduction  //
   
endfunction


function v = to_daily (u)

   ## usage: v = to_daily (u)
   ## 
   ## transform to daily

   if !isfield(u, "id")
      for [vu ku] = u
	  v.(ku) = to_daily (vu) ;
      endfor
      return ;
   endif

   v = u ;
   if !all(diff(datenum(v.id)) == 1) && length(unique(v.id(:,4))) > 1
      if !isempty(strfind(v.name, "precip"))
	 [v.id, v.x] = dayavg(v.id, v.x, 24) ;
      else
	 [v.id, v.x] = dayavg(v.id, v.x, 0) ;
      endif
   endif
   v = fill_date(v) ; # fix broken scenarios

endfunction


function v = to_grid (u, un, iptr)

   ## usage: v = to_grid (u, un, iptr)
   ## 
   ## regrid to common resolution for predictor iptr
   ## UGLY

   global PAR
   tol = 1e-3 ;

   v = u ;

   Jlon = PAR.lon(1) <= v.lon & v.lon <= PAR.lon(2) ;
   Jlat = PAR.lat(1) <= v.lat & v.lat <= PAR.lat(2) ;

   if numel(PAR.nlon) > 0 && numel(PAR.nlat) > 0 && numel(PAR.gausslat) > 0
      [lon lat] = target_ll(PAR.nlon{iptr}, PAR.nlat{iptr}, PAR.gausslat{iptr}) ;
   else
      lon = v.lon(Jlon) ;
      lat = v.lat(Jlat) ;
   endif

   if chklonlat(unique(v.lon(Jlon)), unique(lon), unique(v.lat(Jlat)), unique(lat), tol)
      if v.stdgrd
	 if !isfield(v, "Llat") || v.Llat
	    [wlon wlat] = meshgrid(v.lon, v.lat) ;
	 else
	    [wlat wlon] = meshgrid(v.lat, v.lon) ;
	 endif
	 Jptr = (PAR.lon(1)<=wlon & wlon <= PAR.lon(2) & PAR.lat(1) <= wlat & wlat <= PAR.lat(2)) ;
	 v.lon = v.lon(any(Jptr, 1)) ;
	 v.lat = v.lat(any(Jptr, 2)) ;
      else
	 Jptr = Jlon & Jlat ;
	 v.lon = v.lon(Jlon) ;
	 v.lat = v.lat(Jlat) ;
      endif
      v.x = v.x(:, Jptr) ;
   else
      if v.stdgrd
	 qlon = mean(abs(diff(lon))) / mean(abs(diff(v.lon))) ;
	 qlat = mean(abs(diff(lat))) / mean(abs(diff(v.lat))) ;
	 if ((qlon > 0 && qlon < PAR.qint) || (qlat > 0 && qlat < PAR.qint))
	    warning("xds:xds", "ptr: interpolating %s to finer grid: qlon=%f, qlat=%f (PAR.qint=%f)", un, qlon, qlat, PAR.qint) ;
	    exit(2) ;
	 endif
      endif
      printf("ptr: regridding %s using %dx%d\n", un, length(lon), length(lat)) ; 
      v = regrid(v, lon, lat) ;
      v.lon = lon' ; v.lat = lat' ; v.x = v.g ;
   endif

endfunction


function [lon lat] = target_ll (nlon, nlat, gausslat)

   ## usage: [lon lat] = target_ll (nlon, nlat, gausslat)
   ## 
   ## calculate target grid

   global PAR

   lon = linspace(0, 360, nlon+1)'(1:end-1) ;
   if gausslat
      lat = gauss2lats(nlat)' ; lat = [-lat ; flipud(lat)] ;
   else
      ## d = (180 / nlat) / 2 ;
      d = 0 ;
      lat = linspace(-90+d, 90-d, nlat+1)' ;
   endif

   lon(lon > 180) = lon(lon > 180) - 360 ;

   lon = lon(PAR.lon(1) <= lon & lon <= PAR.lon(2)) ;
   if (j=find(diff(lon) < 0)) ; # wrap longitudes
      Jlon = [j+1:length(lon) 1:j] ;
      lon = lon(Jlon) ;
   endif
   lat = lat(PAR.lat(1) <= lat & lat <= PAR.lat(2)) ;

endfunction


function v = do_norm (v, cfile)

   ## usage: do_norm (v, cfile)
   ## 
   ## calculate climatologies

   global PAR PRJ isCUR CDIR

   file2 = fullfile(CDIR, "ptr.ob") ;
   if isCUR && !isnewer(cfile, file2) # calc. annual cycle
      vw = selper(v, sdate(v.id, PAR.cper(1,:), PAR.cper(2,:))) ;
      if !PAR.anc
	 [z, clim.xm, clim.xs] = anom(vw.x) ;
      else
	 if v.isP
	    [clim.ps, clim.pp] = anncyclog(vw) ;
	 else
	    [clim.ps, vw] = anncyc(vw) ;
	    clim.pp = prbfit(vw.a) ;
	 endif 
      endif
      save(cfile, "clim") ;
      printf("ptr: --> %s\n", cfile) ; 
   endif

   load(cfile)
   printf("ptr: <-- %s\n", cfile) ; 
   if !PAR.anc
      v.z = anom(v.x, clim.xm, clim.xs) ;
   else
      if v.isP
	 v = anncyclog(v, clim.ps, clim.pp) ;
	 v.z = pfy(v.z, clim.pp(1,:)) ;
      else
	 v = anncyc(v, clim.ps) ;
	 v.z = prbfit(v.a, clim.pp) ;
      endif 
      v = rmfield(v, "a") ; 
   endif 

endfunction


function [pc v] = do_eof (u, avar, svar)

   ## usage: [pc v] = do_eof (u, avar, svar)
   ## 
   ## 

   global PAR isANA isCUR ADIR SDIR keof
   persistent nfun = 0 ; nfun++ ;
   cdata = strrep(SDIR, PAR.sim{:}) ;

   v = u ;

   [w latns] = sort(v.lat, "descend") ;

   Judf = all(isnan(v.z), 1) ;
   if !PAR.udf && any(Judf)
      warning("xds:xds", "ptr: found undefined grid points. Consider setting PAR.udf.") ;
   endif
   if PAR.udf
      efile = fullfile(ADIR, "eof_udf", [avar ".ob"]) ;
      gfile = fullfile(ADIR, "clim_udf/G.ob") ;
      pmkdir(fullfile(ADIR, "udf"), fullfile(ADIR, "eof_udf")) ;
      if !exist(ufile = fullfile(ADIR, "udf", [avar ".ot"]), "file")
	 printf("ptr: --> %s\n", ufile) ;
	 save(ufile, "-text", "Judf") ;
      else
	 ## loading udf mask
	 printf("ptr: <-- %s\n", ufile) ;
	 w = load(ufile) ;
	 Judf = Judf | w.Judf ;
	 printf("ptr: --> %s\n", ufile) ;
	 save(ufile, "-text", "Judf") ;
	 if !isCUR && nfun == 1
	    printf("ptr: ###############################################\n") ;
	    warning("xds:xds", "ptr: Using PAR.udf, so consider the need to re-calibrate!") ;
	    printf("ptr: ###############################################\n") ;
	 endif
      endif
      v.z(:,Judf) = nan ;
   else
      efile = fullfile(ADIR, "eof", [avar ".ob"]) ;
      gfile = fullfile(ADIR, "clim/G.ob") ;
   endif

   if isANA
      sdate(v.id(1,:), PAR.cper(1,:), PAR.cper(2,:)) ;
      vw = selper(v, sdate(v.id, PAR.cper(1,:), PAR.cper(2,:))) ;
      eof.G = covadj(vw.z) ;
      eof.E = pca(vw.z, PAR.trunc) ;
      [eof.nr eof.nc] = size(eof.E) ; eof.lon = v.lon ; eof.lat = v.lat ; eof.latns = latns ;
      printf("ptr: --> %s\n", efile) ; 
      save(efile, "eof") ;
   else
      load(efile) ;
      printf("ptr: <-- %s\n", efile) ;
   endif

   if strfind(PAR.ptradj, "sgl")
      printf("ptr: adjusting covariance for %s\n", svar) ;
      mfile = fullfile(CDIR, "clim", [avar ".m.ob"]) ;
      if isCUR
	 vw = selper(v, sdate(v.id, PAR.cper(1,:), PAR.cper(2,:))) ;
	 wm = mean(vw.z) ;
	 save(mfile, "wm")
	 printf("ptr: --> %s\n", mfile) ; 
	 w = anom(v.z, wm) ;
      else
	 if strfind(PAR.ptradj, ".loc")
	    [w wm] = anom(v.z) ;
	 else
	    load(mfile) ;
	    printf("ptr: <-- %s\n", mfile) ; 
	    w = anom(v.z, wm) ;
	 endif
      endif 
      w = covadj(w, eof.G) ;
      v.z = repmat(wm, rows(w), 1) + w ;
   endif

   printf("ptr: EOFs: %s  %d  (%d to %d)\n", svar, eof.nc, keof+1, keof+=eof.nc) ;

   [eof.E lp lmap] = latp(eof.E, eof.latns, latns) ;
   if lmap
      printf("latp: matching latitudes:\n") ;
      printf(["latp: " repmat("%7.2f", 1, length(v.lat)) "\n"], eof.lat(lp)) ;
   endif

   d1 = ndims(v.z) ; d2 = ndims(eof.E) - 1 ;
   pc = mult(v.z, eof.E, d1, d2) ;

endfunction


function res = upd_PC (PC, v, x, iptr)

   ## usage:  res = upd_PC (PC, v, x, iptr)
   ##
   ## update PC from ptr

   if isempty(PC.id)
      res = v ;
      res.x = x ;
   else
      [IPC Iv] = common(PC.id, v.id) ;
      res.id = v.id(Iv,:) ;
      res.x = [PC.x(IPC,:) x(Iv,:)] ;
   endif

   res.J = [PC.J repmat(iptr, 1, columns(x))] ;
   res.cal = v.cal ;

endfunction


function u = mltadj (v)

   ## usage: y = mltadj (v)
   ## 
   ## multivariate covariance adjustment

   global PAR isCUR ADIR CDIR SDIR

   u = v ;

   printf("ptr: adjusting multi-field covariance\n") ;
   mfile = fullfile(CDIR, "pc.m.ob") ;
   if isCUR
      uw = selper(u, sdate(u.id, PAR.cper(1,:), PAR.cper(2,:))) ;
      wm = mean(uw.x) ;
      save(mfile, "wm")
      printf("ptr: --> %s\n", mfile) ; 
      w = anom(u.x, wm) ;
   else
      if strfind(PAR.ptradj, ".loc")
	 [w wm] = anom(u.x) ;
      else
	 load(mfile) ;
	 printf("ptr: <-- %s\n", mfile) ; 
	 w = anom(u.x, wm) ;
      endif
   endif

   printf("ptr: <-- %s\n", gfile = fullfile(ADIR, "G.ob")) ;
   load(gfile) ; Ga = G ;
   printf("ptr: <-- %s\n", gfile = fullfile(CDIR, "G.ob")) ;
   load(gfile) ; Gc = G ;
   w = covadj(w, Ga, Gc) ;
   u.x = repmat(wm, rows(w), 1) + w ;

endfunction


function ptr = read_ptr (gdir)

   ## usage:  ptr = read_ptr (gdir)
   ##
   ## read predictor fields

   global PAR 

   plist = glob({fullfile(gdir, "*.nc"), fullfile(gdir, "*/*.nc")})' ;
   if exist(pfile = fullfile(gdir, "ptr.lst"), "file")
      printf("read_ptr: <-- %s\n", pfile) ;
      plist = union(plist, expand_plist(pfile)) ;
   endif
   if isempty(plist)
      error("read_ptr: no matching files in %s\n", fullfile(gdir, "ptr.lst")) ;
   endif
   [f1 f2 f3] = cellfun(@(f) fileparts(f), plist, "UniformOutput", false) ;
   ## I = eliminate_dups(f2) ; plist = plist(I) ;

   fid = fopen(fullfile(gdir, "done_ptr.lst"), "a+") ;
   ptr = [] ; f = tmpnam ;
   for ncfile = plist

      ncfile = ncfile{:} ;
      printf("read_ptr: <-- %s\n", ncfile) ;

      if isempty(stat(ncfile))
	 nbytes = 0 ; # url
      else
	 nbytes = round(stat(ncfile).size) ; # local
      endif
      if ((k = ceil(nbytes/(0.5*PAR.mem))) > 1)
 	 if system("which ncks", 1)
	    error("ptr: need nco package. Exiting.\n") ;
	 endif 
	 if PAR.dbg, printf("read_ptr: splitting %s into about %d pieces.\n", ncfile, k+1) ; endif
	 ptr = nco_GCM(ncfile, k, ptr) ;
      else
	 w = read_GCM(ncfile) ;
	 if isempty(w), continue ; endif
	 w = to_daily(w) ;
	 ptr = join_ptr(ptr, w) ;
      endif 

      ##printf("read_ptr: --> %s\n", f) ;
      ##save(f, "ptr") ;
      fputs(fid, [ncfile "\n"]) ; fflush(fid) ;

   endfor
   fclose(fid) ;

endfunction


function s = nco_GCM (ncfile, k, r)

   ## usage:  s = nco_GCM (ncfile, k, r)
   ##
   ## read GCM for large files, using nco to split into nbytes

   s = r ;

   [w n] = system(["ncks -m -v time " ncfile " | awk '/time dimension/{print $7}'"], 1) ;
   n = str2num(n) - 1 ;
   n0 = floor(n/k) ;

   i = 0 ;
   while (i <= n)
      tfile = tmpnam ;
      is = num2str(i) ;
      if (i+n0-1 < n)
	 ins = num2str(i+n0-1) ;
      else
	 ins = num2str(n) ;
      endif 
      i += n0 ;
      system(["ncks -Od time," is "," ins " " ncfile " " tfile]) ;
      w = read_GCM(tfile) ;
      # try
	 s = join_ptr(s, w) ;
      # catch
      # 	 exit ;
      # end_try_catch
      system(["rm " tfile]) ;
   endwhile 

endfunction


function s = read_GCM (ncfile)

   ## usage:  s = read_GCM (ncfile)
   ##
   ## read GCM files

   f = dbstack ;

   global PAR

   ncfile = tilde_expand(ncfile) ;

   [nc lDAP] = read_netcdf(ncfile) ;
   if !strcmp(class(nc), "ncfile")
      s = [] ;
      return ;
   endif

   for nv=ncvar(nc)
      nv = nv{:} ;
      if strcmp(ncdatatype(nv), "char"), continue ; endif 
      if (!ncisdim(nc, nv) && !strcmp(nv.long_name, "longitude")
	  && !strcmp(nv.long_name, "latitude"))
	 continue ;
      endif
      ws.(ncname(nv)) = nv(:) ;
   endfor
   STDGRD = false ; j = 0 ;
   for nv=ncdim(nc)
      nv = nv{:} ;
      switch ncname(nv)
	 case {"longitude" "lon"}
	    ws.lon = nc{ncname(nv)}(:) ;
	    ws.nlon = numel(ws.lon) ;
	    STDGRD = true ;
	 case {"latitude" "lat"}
	    ws.lat = nc{ncname(nv)}(:) ;
	    ws.nlat = numel(ws.lat) ;
	    STDGRD = true ;
	 otherwise
	    ws.(["n" ncname(nv)]) = nv(:) ;
      endswitch
   endfor
   if isfield(ws, "lon")
      J = ws.lon > 180 ;
      ws.lon(J) = ws.lon(J) - 360 ;
   endif
   if isfield(ws, "longitude"),  ws = rmfield(ws, "longitude") ; endif
   if isfield(ws, "latitude"),  ws = rmfield(ws, "latitude") ; endif

   jLev = find(index(fieldnames(ws), "lev")) ;
   if hasLev = (!isempty(jLev))
      jLev = jLev(1) ; lev = fieldnames(ws){jLev} ;
      if strcmp(nc{lev}.units, "Pa"), ws.(lev) = ws.(lev) / 100 ; endif 
   endif

   hasTime = false ;
   if isfield(ws, "time")
      wdate = "time" ; wtime = "date" ;
   endif 
   if isfield(ws, "date")
      wdate = "date" ; wtime = "time" ;
      hasTime = isfield(ws, "time") ;
   endif 
   try
      if isempty(wcal = nc{(wtime)}.calendar)
	 wcal = "gregorian" ;
      endif 
   catch
      if isempty(wcal = nc{(wdate)}.calendar)
	 wcal = "gregorian" ;
      endif 
   end_try_catch

   wstr = strsplit(nc{(wdate)}.units, " ", STRIP_EMPTY=true) ;
   if strcmp(wstr(1){:}, "hours")
      wq = 24 ;
   else
      wq = 1 ;
   endif

   clear id
   if strcmp(wstr(2){:}, "since")
      d0 = [sscanf(wstr(3){:}, "%d-%d-%d"); 0]' ;
      if (strcmp(wcal, "gregorian") || strcmp(wcal, "standard")) && date_cmp(d0(1:3), [1582 10 14])
	 jd = nc{(wdate)}(:)/wq - 2 ;
      else
	 jd = nc{(wdate)}(:)/wq ;
      endif
      j0 = date2cal(d0, wcal) ;
      [id(:,1), id(:,2), id(:,3), id(:,4)] = cal2date(j0 + jd, wcal) ;
   else
      id = dat(nc{"time"}(:)) ;
   endif

   lon0 = PAR.lon(1) ; lon1 = PAR.lon(2) ;
   lat0 = PAR.lat(1) ; lat1 = PAR.lat(2) ;

   ## widen selection area for common grids
   d = (lon1 - lon0) * PAR.qgrd ;
   lon0 = lon0 - d ; lon1 = lon1 + d ; lat0 = lat0 - d ; lat1 = lat1 + d ;

   for nv=ncvar(nc)
      nv = nv{:} ;
      if ncisdim(nc, nv), continue ; endif 
      k = ncname(nv) ;
      ## if isempty(nv._FillValue) && isempty(nv.missing_value), continue ; endif # neglect non-fields
      if length(ncdim(nv)) < 3, continue ; endif # neglect non-fields
      if (!STDGRD && (strcmp(k, "lon") || strcmp(k, "lat"))), continue ; endif

      if STDGRD
	 Jx = find(lon0 <= ws.lon & ws.lon <= lon1) ;
	 Jy = find(lat0 <= ws.lat & ws.lat <= lat1) ;
      else
	 [ll0, ll1] = rotfind([lon0  lat0], [lon1 lat1], ws.lon, ws.lat) ;
	 Jx = ll0(1):ll1(1) ;
	 Jy = ll0(2):ll1(2) ;	  
	 if PAR.dbg
	    c0 = min([ws.lon(Jx, Jy)(:), ws.lat(Jx, Jy)(:)]) ;
	    c1 = max([ws.lon(Jx, Jy)(:), ws.lat(Jx, Jy)(:)]) ;
	    printf("ptr: nonstandard grid with range (%4.1f,%4.1f) (%4.1f,%4.1f)\n", c0, c1) ;
	 endif
      endif

      if (j=find(diff(Jx) > 1)) ; # wrap longitudes
	 Jx = Jx([j+1:length(Jx) 1:j]) ;
      endif

      if isempty(Jx) || isempty(Jy)
	 error("xds:xds", "%s: no ptr values, check lon, lat", f(1).name) ;
      endif

      idx = {':'}(ones (1, ndims(nv))) ; 
      if STDGRD
	 wlon = ws.lon(Jx) ; wlat = ws.lat(Jy) ;
	 [ilon, ilat] = inddim(nv) ;
	 idx{ilon} = Jx ; idx{ilat} = Jy ;
      else
	 wlon = ws.lon(Jx,Jy) ; wlat = ws.lat(Jx,Jy) ;
	 idx{end-1} = Jx ; idx{end} = Jy ;
      endif

      if lDAP
	 w = cat_idx(nv, idx) ; # ugliest hack ever (waiting for improved octcdf)
      else
	 w = nv(:) ; N = ones(1, ndims(w)) ; w = w(idx{:}) ;
      endif
      N(1:ndims(w)) = size(w) ;
      w = squeeze(reshape(w, [N(1:end-2) prod(N(end-1:end))])) ;

      # w = nv(idx{:}) ; N = size(w) ; FIXME: wait for logical indexing in octcdf

      try
	 w(w == nv._FillValue) = NaN ;
      end_try_catch 
      try
	 w(w == nv.missing_value) = NaN ;
      end_try_catch 
      try
	 w = nv.scale_factor * w + nv.add_offset ;
      end_try_catch 

      ## expand times, levels...
      dl = 2 ;
      if hasTime, dl += 1 ; endif 
      if hasLev
	 for j=1:ws.(["n" lev])
	    if isfield(PAR, "lev") && !isempty(PAR.lev) && isfield(ws, "lev") && !isempty(ws.lev) && !(abs(ws.lev(j)-PAR.lev)/PAR.lev< 1e-3), continue, endif
	    if isfield(PAR, "plev") && !isempty(PAR.plev) && isfield(ws, "plev") && !isempty(ws.plev) && ws.plev(j) != PAR.plev, continue, endif
	    kk = [k "_" strrep(num2str(ws.(lev)(j)), ".", "_")] ;
	    s.(kk).id = id ;
	    s.(kk).x = nsqueeze(w, dl, j) ;

	    s.(kk).x = outlier(s.(kk).x) ;

	    s.(kk).time = ws.time ;
	    s.(kk).ntime = ws.ntime ;
	    s.(kk).cal = wcal ;
	    s.(kk).lon = wlon ; s.(kk).lat = wlat ;
	    if STDGRD, s.(kk).Llat = ilon > ilat ; endif
	    s.(kk).name = nv.long_name ;
	    s.(kk).stdgrd = STDGRD ;
	    if hasTime, s.(kk).N = N ; endif
	 endfor
      else
	 s.(k).id = id ;
	 s.(k).x = squeeze(w) ;

	 s.(k).x = outlier(s.(k).x) ;

	 s.(k).time = ws.time ;
	 s.(k).ntime = ws.ntime ;
	 s.(k).cal = wcal ;
	 s.(k).lon = wlon ; s.(k).lat = wlat ;
	 if STDGRD, s.(k).Llat = ilon > ilat ; endif
	 s.(k).name = nv.long_name ;
	 s.(k).stdgrd = STDGRD ;
	 if hasTime, s.(k).N = N ; endif
      endif

      if !STDGRD
	 s.(k).id(:,4) = 12 ;
	 # s.(k) = fill_date(s.(k)) ;
      endif

   endfor

   if hasTime
      s = convhasTime(s) ;
   endif

endfunction


function s = convhasTime (r)

   ## usage:  s = convhasTime (r)
   ##
   ## convert times when time dimension

   s = r ;

   for [v k] = r

      if !isfield(v, "id"), continue ; endif 
      if !isfield(v, "time"), continue ; endif 
      if !isfield(v, "ntime"), continue ; endif 

      for i=1:rows(v.id)
	 for j=1:length(v.time)
	    id((i-1)*v.ntime+j,:) = [v.id(i,1:3) v.time(j)] ;
	 endfor
      endfor
      v.id = id ;

      nd = length(v.N) ;
      if (nd > 2)
	 v.x = permute(v.x, [3:nd, 2, 1]) ;
	 N = size(v.x) ;
	 v.x = reshape(v.x, [N(1:nd-2), N(nd-1)*N(nd)]) ;
	 v.x = permute(v.x, [nd-1,1:nd-2]) ;
      endif 
      v = rmfield(v, "time") ;
      v = rmfield(v, "ntime") ;
      s.(k) = v ;

   endfor 

endfunction


function res = ncisdim (nc, nv)

   ## usage:  res = ncisdim (nv)
   ##
   ## netcdf variable nv of nc is dim

   res = false ;
   for nd = ncdim(nc)
      res = res || strcmp(ncname(nd{:}), ncname(nv)) ;
   endfor 

endfunction


function id = dat(idat, j0)

%usage:  id = dat(idat)

   iw = idat;
   if (nargin < 2), j0 = 4; endif 
   ind=ones(rows(idat),j0);
   id=ones(rows(idat),3);

   for j=1:j0
      n = 10^(2*(j0-1-j));
      ind(:,j) = iw(:)./n;
      iw = rem(iw,n);
      id(:,j) = fix(ind(:,j)) ;
   endfor

   if (j0 == 4)
      id(:,j0) = id(:,j0)/100 * 24 ;
   endif 

endfunction


function [ilon, ilat] = inddim (nv)

   ## usage:  [ilon, ilat] = inddim (nv)
   ##
   ## find longitudes and latitudes

   i = 0 ;
   for d = ncdim(nv)
      i++ ;
      if strfind(ncname(ncdim(nv){i}), "lon"), ilon = i ; endif
      if strfind(ncname(ncdim(nv){i}), "lat"), ilat = i ; endif
   endfor
   
endfunction


function res = cat_idx (nv, idx)

   ## usage:  res = cat_idx (nv, idx)
   ##
   ## concatenate netcdf arrays

   for j=1:numel(idx)

      if strcmp(J = idx{j}, ":") 
   	 widx{j}{1} = J ;
   	 continue
      endif ;
      kJ = find(diff(J) < 0) ;

      if isempty(kJ)
   	 widx{j}{1} = [num2str(J(1)) ":" num2str(J(end))] ;
      else
   	 widx{j}{1} = [num2str(J(1)) ":" num2str(J(kJ))] ;
   	 widx{j}{2} = [num2str(J(kJ+1)) ":" num2str(J(end))] ;
      endif

   endfor

   for j = 1:numel(idx)
      n(j) = numel(widx{j}) ;
   endfor

   res = [] ; cI = cell(numel(idx),1) ; I0 = zeros(1, numel(idx)) ;
   for j = 1:prod(n)
      [cI{:}] = ind2sub(n, j) ; I = cell2mat(cI)' ;
      J = "" ;
      for i = 1:numel(idx)
   	 J = [J widx{i}{I(i)} ","] ;
      endfor
      J = J(1:end-1) ;
      eval(["w = nv(" J ") ;"]) ;

      i0 = find(I != I0)(1) ; I0 = I ;
      res = cat(i0, res, w) ;
   endfor
 
endfunction


function [nc lDAP] = read_netcdf (ncfile)

   ## usage:  [nc lDAP] = read_netcdf (ncfile)
   ##
   ## read netcdf file, from OPeNDAP or fileserver

   pkg load octcdf netcdf
   if exist(cert = fullfile(pwd, "cert.sh"))
      system(["sh " cert]) ;
   endif

   try
      ncinfo(ncfile) ;
   catch
      warning("xds:xds", "could not read: %s", ncfile) ;
   end_try_catch

   f = dbstack ;
   global PAR

   if isempty(strfind(ncfile, "://"))

      lDAP = false ;
      if stat(ncfile).size == 0
	 warning("xds:xds", "%s: empty file: %s", f(1).name, ncfile) ;
	 unlink(ncfile) ;
	 nc = [] ;
      else
	 nc = netcdf(ncfile, "r") ;
      endif

   else

      ##try
	 lDAP = true ;
	 ##f = tmpnam ;
	 ##wget_args = [" --certificate=" PAR.credsdir "/credentials.pem --private-key=" PAR.credsdir "/credentials.pem \
	 ##--ca-directory=" PAR.credsdir "/certificates -qO " f] ;
	 ##[st out] = system(cmd = ["wget " wget_args " " [ncfile ".dds"]]) ;
	 ##unlink(f) ;
	 for j=1:10		# try hard
	    [st out] = system(["ncdump -h " ncfile]) ;
	    if st == 0
	       nc = netcdf(ncfile, "r") ;
	       break ;
	    else
	       warning("xds:xds", "read_netcdf: trouble reading %s", ncfile) ;
	       if exist(cert = fullfile(pwd, "cert.sh"))
		  system(["sh " cert]) ;
	       endif
	       sleep(5) ;
	    endif
	 endfor
	 if st ~= 0
	    warning("xds:xds", "read_netcdf: not read: %s", ncfile) ;
	    lDAP = false ; nc = NaN ;
	    return ;
	 endif
      ##catch
      ##	 lDAP = false ;
      ##	 [f1 f2 f3] = fileparts(ncfile) ;
      ##	 [info, err, msg] = stat(f = fullfile("data", PAR.sim{1}, [f2 f3])) ;
      ##	 if err
      ##	    wget_args = [" --certificate=" PAR.credsdir "/credentials.pem --private-key=" PAR.credsdir "/credentials.pem \
      ##	 --ca-directory=" PAR.credsdir "/certificates -qO " f] ;
      ##	    ncinfo(ncfile)
      ##	    [err out] = system(cmd = ["wget " wget_args " " ncfile]) ;
      ##	 endif
      ##	 if err
      ##	    error("ptr:read_netcdf: file not found: %s\nexecuting: %s\nstatus: %d: %s", ncfile, cmd, st, out) ;
      ##	 else
      ##	    nc = netcdf(f, "r") ;
      ##	 endif
      ##end_try_catch

   endif

endfunction


function [r0, r1] = rotfind (ll0, ll1, lon, lat)

   ## usage: [r0, r1] = rotfind (ll0, ll1, lon, lat)
   ## 
   ## find rotated grid section

   I = lon >= ll0(1) & lon <= ll1(1) & lat >= ll0(2) & lat <= ll1(2) ;
   [i j] = find(I) ;
   [r0(1) r0(2)] = deal(min(i), min(j)) ;
   [r1(1) r1(2)] = deal(max(i), max(j)) ;
	 
endfunction


function v = regrid (u, lon, lat)

   ## usage:  v = regrid (u, lon, lat)
   ##
   ## regrid u.x with u.lon, u.lat, using vectors lon, lat

   global PAR

   v = u ;

   nr = rows(u.x) ;

   if v.stdgrd
      nlat = length(u.lat) ; nlon = length(u.lon) ;
      u.x = reshape(u.x, nr, nlat, nlon) ;
      [u.lon plon] = sort(u.lon) ;
      [u.lat plat] = sort(u.lat) ;
      u.x = u.x(:, plat, plon) ;
      c = mat2cell(u.x, ones(1, nr), nlat, nlon) ;
   else
       [nx ny] = size(u.lon) ; u.x = reshape(u.x, [], nx, ny) ;
       c = mat2cell(u.x, ones(1, nr), nx, ny) ;
   endif

   y = parcellfun_(@(x) interp2(u.lon, u.lat, squeeze(x), lon, lat'), c, PAR.par_opt{:}) ;
   y = cell2mat(y) ;
   y = reshape(y, length(lat), nr, length(lon)) ;
   y = permute(y, [2 1 3]) ;
   v.g = reshape(y, nr, length(lat)*length(lon)) ;
   
endfunction


function v = common_def (u)

   ## usage: v = common_def (u)
   ## 
   ## 

   v = u ;

   I = ID = cell(1, numfields(u)) ;
   j = 0 ;
   for [s k] = u
       ID{++j} = s.id ;
   endfor

   [I{:}] = common(ID{:}) ;

   j = 0 ;
   for [s k] = u
       j++ ;
       if isfield(s, "id"), v.(k).id = s.id(I{j},:) ; endif
       if isfield(s, "x"), v.(k).x = s.x(I{j},:) ; endif
       if isfield(s, "z"), v.(k).z = s.z(I{j},:) ; endif
       if isfield(s, "a"), v.(k).a = s.a(I{j},:) ; endif
   endfor

endfunction


function L = chklonlat (lon1, lon2, lat1, lat2, tol)

   ## usage: L = chklonlat (lon1, lon2, lat1, lat2, tol)
   ## 
   ## check lon lat for equality

   L = numel(lon1) == numel(lon2) && numel(lat1) == numel(lat2) ;

   if !L, return, endif

   L = norm(lon1 - lon2) <= tol && norm(lat1 - lat2) <= tol ;

endfunction
